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ABSTRACT 


This  thesis  will  examine  the  effect  of  nonlinearities  on  the  propagation  of  orbit 
uncertainties  in  order  to  gain  insight  into  the  accurateness  of  the  estimation  of  covariance 
with  time.  Many  real-world  applications  rely  on  a  first-order  approximation  of  nonlinear 
equations  of  motion  for  propagation  of  orbit  uncertainty.  The  nonlinear  effects  that  are 
ignored  during  the  linearization  process  can  greatly  influence  the  accuracy  of  the 
solution.  A  comparative  analysis  of  linear  and  nonlinear  orbit  uncertainty  propagation  is 
presented  in  order  to  attempt  to  detennine  when  linearized  uncertainty  becomes  non- 
Gaussian.  An  examination  of  performance  metrics  is  then  accomplished  to  compare 
linearly  propagated  uncertainty  to  uncertainty  propagated  using  a  second-order 
approximation.  An  attempt  is  then  made  to  develop  a  perfonnance  metric  that  determines 
when  propagated  uncertainty  is  no  longer  Gaussian.  The  results  show  it  is  difficult  to 
determine  a  clear  method  of  defining  when  the  linear  approximated  uncertainty  is  no 
longer  Gaussian,  but  there  are  metrics  that  can  be  implemented  given  a  user-defined 
threshold  of  performance. 
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I.  INTRODUCTION 


This  thesis  will  examine  the  effect  of  nonlinearities  on  the  propagation  of  orbit 
uncertainties  in  order  to  gain  insight  into  the  accurateness  of  covariance  estimation  with 
time.  Accurate  covariance  can  be  a  valuable  tool  in  space  situational  awareness. 
Covariance  infonnation  can  be  used  in  the  correlation  of  uncorrelated  tracks  (UCTs) 
[  1  ]— [3],  computing  probability  of  collision  [4]  and  sensor  tasking  [5]. 

Air  Force  Space  Command  (AFSPC)  currently  tracks  over  22,000  orbiting  objects 
via  the  Air  Force  Space  Surveillance  Network  (SSN).  When  the  SSN  detects  an  object 
that  does  not  correlate  to  an  object  in  the  Space  Object  Catalog,  it  is  considered  a  UCT. 
Dr.  Alfriend  proves  in  [1]  that  the  use  of  covariance  for  uncorrelated  track  association  is 
the  optimal  measure  as  long  as  the  uncertainties  remain  Gaussian.  Provided  the 
covariance  is  accurate,  it  can  then  be  used  to  associate  the  UCT  in  a  more  statistically 
valid  and  automated  manner  than  the  current  labor-intensive  process  [2],  In  addition  to 
UCT  association,  covariance  is  also  useful  in  calculating  the  probability  of  collision 
between  orbiting  objects. 

The  increased  reliance  on  space  systems  in  recent  history  has  led  to  increasing 
probability  of  collisions  of  on-orbit  assets.  A  less  accurate  method  of  the  collision 
avoidance  calculation  involves  predicting  conjunctions  that  fall  within  boxes  centered  on 
the  two  orbiting  objects  in  question.  The  disadvantage  of  this  is  that  it  does  not  take  into 
consideration  the  accuracy  of  the  predicted  boxes.  Prior  to  launch  of  the  International 
Space  Station,  NASA  began  using  covariance  for  the  probability  of  collision  calculations. 
Other  Department  of  Defense  and  national  agencies  have  begun  to  use  the  same  approach 
for  predicting  conjunctions  and  the  ensuing  probability  of  collision  of  on-orbit  assets. 
Much  like  UCT  association,  an  accurate  covariance  is  a  valuable  tool  in  determining  with 
greater  precision  the  probability  of  collision  between  two  objects  that  will  pass  close  to 
each  other  [4].  In  tenns  of  determining  whether  spacecraft  collision  avoidance 
maneuvers  are  necessary,  the  accuracy  of  the  probability  of  collision  calculation  can  then 
be  a  key  factor  with  respect  to  propellant  conservation  and  on-orbit  safety. 
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Most  applications  involving  uncertainty  propagation  involve  a  linear 
approximation  of  a  mathematical  model  of  the  nonlinear  equations  of  motion.  The 
linearization  of  the  equations  of  motion  for  uncertainty  propagation  simplifies  things  a 
great  deal;  however  the  nonlinear  effects  that  are  ignored  during  this  linearization  process 
can  greatly  influence  the  accurateness  of  the  solution.  With  this  in  mind,  this  thesis 
primarily  attempts  to  answer  two  main  questions.  First,  can  examining  the  perfonnance 
of  nonlinear  uncertainty  propagation  provide  any  insight  into  when  the  linear 
approximation  fails?  Second,  which  errors  in  the  initial  orbit  uncertainty  cause  the 
propagated  covariance  to  be  non-Gaussian  and  how  large  do  these  errors  have  to  be? 

Through  comparative  analysis  of  linear  and  nonlinear  orbit  uncertainty 
propagation  techniques,  we  can  attempt  to  determine  when  the  effects  of  the 
nonlinearities  cause  the  inaccuracy  of  the  covariance  to  be  large  enough  to  render  it  no 
longer  useful.  In  order  to  gain  insight  into  when  the  linear  approximation  fails,  this  thesis 
will  consider  the  classical  orbital  elements  and  two-body  motion  with  a  comparison  of  a 
first-order  approximation  to  a  second-order  approximation. 

Chapter  II  is  a  background  chapter.  It  gives  a  brief  review  of  the  Clohessy- 
Wiltshire  equations,  the  definition  of  the  orbit  state  that  will  be  used  and  then  an 
overview  of  the  higher  order  approximation  that  will  be  used  in  the  comparative  analysis. 
Chapter  III  presents  the  results  of  the  comparative  analysis  and  Chapter  IV  addresses  the 
topic  of  determining  appropriate  metrics  to  assess  the  performance  of  a  first-order 
approximation. 
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II.  BACKGROUND 


A.  FORMULATION  OF  THE  PROBLEM 

1.  Overview 

Expressions  that  describe  the  dynamics  of  orbital  motion  are  based  on  a  set  of 
nonlinear  differential  equations.  However,  in  many  real-world  applications,  the 
propagation  of  uncertainty  is  accomplished  using  a  first-order  approximation.  Given 
initial  Gaussian  uncertainty,  the  propagated  uncertainty  remains  Gaussian  as  long  as  the 
system  remains  linear.  The  system  however  is  nonlinear,  and  eventually  the  neglected 
nonlinearities  influence  the  accuracy  of  the  propagated  covariance. 

The  dynamics  of  uncertainty  propagation  are  often  described  using  a  first-order 
approximation.  More  complex  models  exist;  however  selecting  the  appropriate  model  is 
a  balance  between  simplicity  and  accuracy.  With  increased  accuracy,  comes  higher  cost 
in  terms  of  computation  time,  computing  resources  and  overall  problem  complexity. 
Insight  into  when  a  particular  mathematical  model  is  no  longer  valid  or  when  it  no  longer 
meets  requirements  is  very  useful.  Determining  this  threshold  will  be  the  primary  focus 
of  this  thesis.  Previous  work  by  Fujimoto  et  al.  [6],  and  Park  and  Scheeres  [7],  has 
shown  significant  increases  in  model  accuracy  can  be  obtained  by  using  a  second-order 
approximation.  Much  of  the  comparative  analysis  presented  in  this  thesis  relies  on  a 
comparison  of  second-order  to  first-order  approximation  for  covariance  propagation. 

The  background  information  presented  in  the  remainder  of  this  chapter  consists  of 
a  brief  review  of  the  Clohessy-Wiltshire  (CW)  equations,  the  definition  of  the  orbit  state, 
and  then  presents  an  overview  of  the  higher  order  approximation  that  will  be  used  for 
comparative  analysis. 

2.  Linearized  Model 

A  linearized  mathematical  model  is  useful  for  describing  the  motion  of  an  orbiting 
vehicle  relative  to  a  reference  orbit.  The  well-known  Clohessy-Wiltshire  equations 
provide  a  linearized  model  that  approximates  relative  orbital  motion.  The  CW  equations 
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are  the  linearization  of  the  equations  of  motion  about  a  circular  orbit  and  spherical  Earth. 
Although  there  are  more  complicated  linearized  models  that  consider  non-circular  orbits 
and  perturbations  [8],  the  CW  equations  can  be  useful  for  demonstrating  the  types  of 
relative  motion  obtained  by  closely  orbiting  objects.  The  derivations  of  the  CW 
equations  can  be  found  in  found  in  [9].  For  purposes  of  this  thesis,  the  CW  equations  are 
merely  stated.  The  CW  equations  are: 

x-2ny-3rrx  =  0 

y  +  2nx  =  Q  (1-1) 

z  +  n2z  =  0 

where  n  is  the  mean  motion  defined  by: 


n  - 


(1.2) 


The  solutions  for  the  CW  equations  are  given  by  the  following: 


V  n  J 


x  =  (4-3cos(«t))x0  +  —  |sin(«t)  +  [  — |(l-cos(«t)) 


y  =  6  (sin  (nt)  -  nt}  x0  +  2 


(x  \ 
Ao 

V  n 


cos 


(. nt)  +  y0+  ^-J(4sin(«t)-3«t)  (1.3) 


z  =  z0cos(nt)  +  |  —  Isin(nt) 
n 


The  CW  frame  is  illustrated  in  Figure  1.  The  x-axis  or  radial  component  lies 
along  the  object’s  radius  vector.  The  y-axis  or  in-track  component  lies  in  the  orbit  plane 
perpendicular  to  the  x-axis.  The  z-axis  or  cross-track  component  is  perpendicular  to  the 
orbit  plane.  The  origin  of  the  CW  frame  moves  with  the  reference  object  and  is  a  rotating 
frame.  The  angular  velocity  of  the  rotating  frame  is  equal  to  the  orbital  angular  velocity. 
This  coordinate  system  is  also  referred  to  as  the  radial,  in-track,  cross-track  or  RIC  frame 
and  is  very  useful  when  discussing  differences  between  two  orbits. 
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The  motion  of  an  anecdotal  low  earth  orbit  object  under  the  CW  model  is 
illustrated  in  Figure  2.  One  can  see  that  the  motion  in  the  radial  and  cross-track 
directions  are  periodic  and  do  not  grow  in  magnitude  under  the  linear  assumption.  The 
in-track  motion  does  change  significantly  with  time.  It  is  the  growth  in  the  in-track 
direction  that  has  higher  order  effects  on  the  linear  model  and  causes  the  uncertainty  to 
no  longer  be  Gaussian.  The  goal  is  to  determine  when  this  happens  as  a  function  of  time 
and  initial  uncertainties.  The  secular  term  from  Equation  (1.3)  is  given  by: 

Secular  Motion  =  (-6nx0  -3  y0)t  (1-4) 


In  terms  of  differential  motion  in  the  orbital  elements  in  the  CW  model, 
inclination,  right  ascension  of  the  ascending  node  and  eccentricity  all  cause  periodic 
variations.  The  in-track  secular  motion  is  caused  by  differentials  in  the  mean  anomaly 
and  semi-major  axis.  In  an  attempt  to  isolate  the  non-Gaussian  behavior  caused  by  the 
secular  motion  in  the  in-track  direction,  the  two-dimensional  state  of  mean  anomaly  and 
semi-major  axis  will  be  used.  This  2D  state  should  capture  the  essential  elements  of  the 
nonlinear  behavior  that  causes  the  non-Gaussian  behavior  of  the  uncertainty. 
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Motion  described  by  Clohessy-Wiltshire  Equations 


x  10*  x  108 


0  rbits 


Figure  2.  Motion  Described  by  CW  Equations 


3.  Definition  of  the  State 

In  an  attempt  to  isolate  the  non-Gaussian  behavior  caused  by  the  secular  motion 
in  the  in- track  direction,  the  state  considered  for  this  thesis  will  consist  of  the  semi-major 
axis  and  mean  anomaly.  The  dynamics  of  an  n  dimensional  system  can  be  expressed  as: 

(1.5) 

xvx2,...,xn) 


where  the  solution  to  Equation  (1.5)  is  given  by: 

x(t)  =  ^(t,x0,t0)  (1.6) 


and  x0  is  the  initial  state  at  time  t0.  Given  the  reference  solution  x (/,/„),  the  state 
deviation  vector  can  be  represented  by: 

£x  =  x(7)-x(r)  (1.7) 


It  follows  then  that: 

Si  =  f(t,x)-f(t,i)  (1.8) 


The  analysis  presented  will  consist  of  the  two  dimensional  state  consisting  of  the 
change  in  semi-major  axis  and  mean  anomaly.  The  “truth”  state  consists  of  two-body 
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motion  with  no  perturbations.  The  equations  of  motion  for  the  two  dimensional  state  are 
expressed  in  Equation  (1.9): 


a  =a0 
a  =  0 


(1.9) 


The  highest  order  approximation  used  for  this  analysis  will  be  to  second-order. 
Although  a  higher  order  approximation  does  provide  a  more  accurate  solution,  the  most 
significant  perfonnance  gain  over  the  linear  approximations  comes  from  the  second-order 
approximation.  To  obtain  the  equations  that  approximate  the  motion,  we  will  expand  in  a 
Taylor  series  about  the  reference  orbit  to  second-order. 


The  state  is  given  by: 


8x  = 


8a 

8M 


(1.10) 


where  the  notation  indicates  a  normalized  dimensionless  term.  The  non-dimensional 
change  in  semi-major  axis  is  given  by: 


a 


(1.11) 


where  the  “bar”  notation  denotes  the  value  on  the  reference  orbit.  The  time  scale  used 
will  be  expressed  in  terms  of  number  of  orbits,  r .  The  time  scale  is  given  by: 


T 


nt 
2  n 


(1.12) 


The  expressions  for  the  change  in  semi-major  axis  and  mean  anomaly  (expanded  to 
second-order)  become: 
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r>  *  c*  * 

da  =  da n 


8M  = 


~  c*  1  ~  *  \- 

-inoa  -\ - ^c?a  j 


r  +  SMn 


(1-13) 


Finally,  using  the  first-order  term  from  Equation  (1.13),  the  linearized  state  transition 
matrix  is: 


®(r) 


'  1  0^ 
k-3ttt  1  y 


(1.14) 


4.  Nonlinear  Uncertainty  Mapping 


Equation  (1.14)  is  the  familiar  linear  state  transition  matrix  determined  by  a  first- 
order  approximation.  In  order  to  capture  the  nonlinear  behavior  neglected  in  the  linear 
approximation,  Park  and  Scheeres  [10]  proposed  a  higher  order  method  of  approximating 
an  orbit  state  and  covariance  by  incorporating  the  higher  order  Taylor  series  terms.  The 
higher  order  state  solutions  as  a  function  of  initial  conditions  are  referred  to  as  state 
transition  tensors  (STT). 


From  Park  and  Scheeres  [10],  the  state  deviation  vector  S\,  the  mean  deviation 
vector  O'm ,  as  well  as  the  uncertainty  P,  expanded  to  second-order  are  given  by: 


=  </>.  d>.  P 

r  i,ar  j^a  a 


0 

aa 


Sm<  =\(<!>i,abPab) 

-  8midmj  +l-</>,JhaP  (PX  +  Pa°X  +  CPL ) 


(1.15) 

(1.16) 
(1.17) 


Where  d'x0  and  P°  are  the  initial  state  and  covariance,  respectively.  Einstein’s 
summation  notation  is  used  here  and  the  subscripts  and  “k”,  respectively  represent  the 
/th  and  kth  component.  Given  the  state  vector  from  Equation  (1.10)  and  (1.13),  the  state 
transition  tensors  expressed  non-dimensionally  are: 


8 


®uM= 

®,«(0= 


dSx 

dSx 

d2Sx 


y-3/TT  1  j 


SSxdSa 
d2Sx 


f  0 
15 

v 


0 


m  0 
2  j 


dSxdSM 


0  0 

vo  oy 


(1.18) 


The  first  of  Equation  (1.18)  represents  the  familiar  linear  state  transition  matrix. 
The  second  and  third  of  Equation  (1.18)  are  the  STTs  and  represent  the  contributions 
from  the  second-order  Taylor  series  expansion.  One  can  see  from  the  expressions  for  the 
second-order  STTs  in  Equation  (1.18),  all  the  higher  order  terms  are  zero  with  the 
exception  of  the  T»,n  term,  which  represents  the  second  partial  derivative  of  SM  with 
respect  to  (Sa  )  . 


<f> 


2,11 


d2SM 
dSa 2 


(1.19) 


Due  to  many  of  the  higher  order  terms  in  the  STTs  being  equal  to  zero,  the 
expressions  for  the  state,  mean,  and  particularly  the  covariance  can  be  simplified  a  great 
deal.  The  expression  for  the  state  given  in  Equation  (1.15)  expanded  from  Einstein’s 
summation  notation  is  written  as: 


5x  = 


(j>x  J  <5xx  +  (j)x  2  Sx°  +  ^(^1  n^Xl^Xl  +  A  12^X\^X2  +  A  2l^X2^Xl  +  A  -,25x\8x^  ) 

^  j  Sxx  +(/)-,  25x°  +— ^2Ii^xi°^xi°  +0,  l2SxxSx2  +<j>2  2  xSx°Sxx  +t/>1  22Sx2Sx2j 


(1.20) 


Given  initial  c)'x0: 


£x0  = 


Sa0 

SMn 


(1.21) 
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when  one  simplifies  Equation  (1.20)  and  discards  the  terms  that  equal  zero,  the  equation 
for  the  state  deviation  vector  as  a  function  of  non-dimensional  time  z  is: 


£x,(r) 


( 

-2>nSa\  + 

V 


Sal 


z  +  SM0 


(1.22) 


Note  that  the  second  term  of  Equation  (1.22)  is  simply  the  second-order  Taylor 
series  expansion  of  SM  shown  in  Equation  (1.13).  In  fact  Equation  (1.22)  and  Equation 
(1.13)  are  equivalent  expressions. 

The  second-order  mean  deviation  vector  from  Equation  (1.16)  is  expressed  as: 


cfrn  = 


+  01,12^1,2  +  01,21-^2,1  +  01,22-^2,2  ) 
+  0242^1,2  +  02,21-^2,1  +  02,22^2,2  ) 


(1.23) 


however,  once  simplified  by  combining  with  Equation  (1.18),  the  mean  shown  in 
Equation  ( 1 .23)  reduces  to: 


(5m(r) 


\5n 


P°T 
1  i,r 


(1.24) 


The  expansion  from  Einstein’s  notation  of  the  second-order  covariance  is  left  for 
the  Appendix.  After  combining  the  expression  for  the  covariance  in  Equation  (1.17)  and 
the  STTs  from  Equation  (1.18),  and  simplifying  terms  equal  to  zero,  the  approximation 
for  the  second-order  covariance  can  be  expressed  as: 


P  = 

y 


KAA 


P° 
"2,1  11 


02,102,1^1  +02,202, 


P° 

2*  22 


Sm2Sm , 


1 

+  — 
4 


>2,„02  ,n(PX 


+  P°P° 
^  ^11^11 


+p; 


>°  p°  \ 
n^ii ) 


(1.25) 


or  as  a  function  of  time  z , 
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P(r) 


P° 

-'ll 

-2>7ttP^ 


-2>7ttP^ 


(1.26) 


Using  Equation  (1.22)  and  (1.26)  we  can  now  propagate  the  state  and  covariance 
using  a  second-order  approximation. 
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III.  COVARIANCE  ANALYSIS 


A.  INTRODUCTION 

This  section  examines  covariance  propagation  through  comparison  of  uncertainty 
ellipsoid  propagation  and  Monte  Carlo  simulations.  The  cumulative  probability 
distribution  for  first-order  and  second-order  approximations  compared  to  the  theoretical 
distribution  for  a  linear  system  is  also  presented.  To  begin,  we  will  examine  how 
uncertainty  propagates  under  the  first-order  and  second-order  approximations. 

B.  ELLIPSE  PROPAGATION 


For  the  uncertainty  ellipse  propagation,  a  demonstration  of  how  an  initial 
Gaussian  distribution  propagates  under  unperturbed  two-body  equations  of  motion  is 
presented.  Only  the  two  dimensional  state  consisting  of  8  a*  and  8M  is  considered. 

For  the  error  ellipsoid  analysis  presented  in  this  section,  106  points  on  the  initial 
3 a  uncertainty  ellipse  are  propagated  using  the  linear  approximation  given  by  the  state 
transition  matrix  in  Equation  (1.14).  The  points  are  also  propagated  using  the  second- 
order  approximation  given  in  Equation  (1.13).  A  Monte  Carlo  simulation  is  then 
performed  by  propagating  106  normally  distributed  points  obtained  from  the  initial 
uncertainty.  The  Monte  Carlo  points  are  propagated  using  the  complete  system 
dynamics: 


C  *  c*  * 

da  =oan 


( 


SM  =  M  -M  = 


A 


(a  +  8a ) 


3t  +  M0 


+  M0 
a 


(2.1) 


Note  that  Equation  (2.1)  is  not  expressed  in  dimensionless  terms  or  in  the 
timescale  z .  The  following  shows  the  expressions  for  the  state  in  non-dimensional  tenns. 
The  non-dimensional  semi-major  axis  is: 
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*  Cl 

a  =  — 
a 


where  the  “bar”  over  the  variable  represents  a  value  on  the  reference  orbit, 
dimensional  time  is  given  by: 


(2.2) 


The  non- 


r  =  ■ 


nt 
2  n 


(2.3) 


Given  that  the  semi-major  axis  is: 


a  =  a  +  da 


(2.4) 


the  normalized  semi-major  axis  is 


a  =l  +  Sa 


(2.5) 


The  expression  for  mean  anomaly  is: 


M  -nt  +  Mn 


(2.6) 


Substituting  t  into  Equation  (2.6)  for  non-dimensional  terms  yields: 


M  =  2tt 


AA  a3 

—  t  +  M0  =  2n.  /— r -r  +  M0 
\  n  )  \l  <r 

3/2 


a  J  t  +  M{ 


(2.7) 


Note  that: 


dM 

dv 


dr 

=  2x(aJm 


(2.8) 


Therefore,  given  that  5M 


M  -M  the  exact  solution  for  the  change  in  mean  anomaly  is: 


SMmt=2x 


(a*yn- 1 


r  +  SMn 


(2.9) 
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1.  Ellipse  Propagation  Results 


By  propagating  the  3cr  uncertainty  ellipsoid  using  first  and  second-order 
propagations  and  comparing  to  the  exact  propagation  of  the  Monte  Carlo  points,  we  are 
able  to  examine  the  validity  of  the  approximations  at  various  propagation  times.  For  this 
scenario,  the  uncertainty  is  given  by: 


(2.10) 


where  the  semi-major  axis  uncertainty  is  21  km,  the  mean  anomaly  uncertainty  is  0.01° 
and  the  reference  semi-major  axis  is  a  =  7000  km.  Given  these  parameters,  the  orbital 
period  is  1.619  hours.  The  standard  deviations  are  given  by: 


The  initial  distribution  as  well  as  the  3cr  ellipse  based  on  the  initial  uncertainty  P0 
are  shown  in  the  top  left  of  Figure  3.  The  remaining  plots  are  the  3cr  ellipse  propagated 
using  first  and  second-order  approximations,  as  well  as  the  Monte  Carlo  points 
propagated  using  the  complete  system  dynamics  at  2,  5,  and  10  orbits. 

The  green  ellipse  represents  the  3cr  ellipse  propagated  using  a  first-order 
approximation  from  Equation  (1.14)  and  the  blue  dashed  ellipse  is  the  second-order 
approximation  from  Equation  (1.13).  The  8a  -  SM  plane  is  rotated  so  the  principal  axes 
of  the  linearly  propagated  ellipse  are  aligned  with  the  horizontal  and  vertical  axis  of  the 
plot.  Ideally,  the  propagated  ellipse  should  continue  to  capture  a  large  number  of  the 
Monte  Carlo  sample  points.  One  can  see  that  the  accuracy  of  the  linearly  propagated 
ellipse  capturing  the  sample  points  significantly  degrades  with  time. 


15 


x  -j o-4 Initial  Error  Distribution  and  3a  Ellipse 


3a  Covariance  After  2  Orbits 


Figure  3.  Initial  Gaussian  Distribution  and  Propagated  Gaussian  Distribution  with 

Corresponding  3<r  Ellipsoids 


The  theoretical  distribution  for  a  two-dimensional  system  is  given  by: 

- k 2 

F{k)  =  \-eY  (2.12) 

where  k  represents  the  number  of  standard  deviations.  For  k=  3,  the  theoretical 
probability  distribution  is  0.9889,  i.e.,  the  3cr  ellipse  should  capture  98.89%  of  the  Monte 
Carlo  sample  points.  As  seen  in  Table  1,  for  the  initial  uncertainty  given  in  (2.1 1),  by  10 
orbits  the  linear  propagation  has  failed  to  represent  a  significant  portion  (nearly  hall)  of 
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the  system  dynamics.  The  illustrations  in  Figure  3  show  the  second-order  propagation 
does  a  much  better  job  of  characterizing  the  full  system  dynamics.  The  results  of  the 
probability  distribution  captured  by  the  first-order  3cr  ellipse  are  summarized  in  Table  1 . 


Theoretical 

Initial 

2  Orbits 

5  Orbits 

10  Orbits 

0.9889 

0.9992 

0.8383 

0.6512 

0.5045 

Table  1.  Linear  Probability  Distribution  through  10  Orbits 


3a  Covariance  After  0.375  Orbits 
0.04  [  •  ’  - 


-10  12  3 

y5M  X  10'4 


3a  Covariance  After  0.75  Orbits 


Figure  4.  Propagated  Gaussian  Distributions  with  Corresponding  3cr  Ellipsoids 
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Although  the  linear  propagation  does  not  fully  represent  the  full  system  dynamics, 
for  a  short  propagation  time,  the  first-order  approximation  may  be  acceptable.  One  can 
see  that  at  2  orbits,  or  3.24  hours,  the  first-order  approximation  still  encompasses  83.8  % 
of  the  sample  points. 

Figure  4  illustrates  the  results  of  various  propagation  times  out  to  1.5  orbits. 
Again,  the  propagated  covariance  is  rotated  so  the  primary  axes  of  the  linear  ellipse 
correspond  to  the  axes  of  the  plot.  One  can  see  via  the  distortion  of  the  second-order 
ellipsoid  and  the  Monte  Carlo  distribution,  there  are  higher  order  effects  even  at  0.375 
orbits,  or  just  over  30  minutes  of  propagation  time.  Despite  some  loss  from  nonlinear 
effects,  depending  on  the  accuracy  required,  the  linear  approximation  at  shorter 
propagation  times  may  be  sufficient.  Table  2  summarizes  the  numerical  results  for  the 
probability  distribution  captured  by  the  linear  approximation. 


Theoretical 

Initial 

0.375  Orbits 

0.75  Orbits 

1.125  Orbits 

1.5  Orbits 

0.9889 

0.9992 

0.9783 

0.9498 

0.9167 

0.8825 

Table  2.  Linear  Probability  Distribution  through  1.5  Orbits 


2.  Ellipse  Rotation 

The  ellipse  rotation  consists  of  rotating  the  major  and  minor  axis  of  the  first-order 
ellipse  to  correspond  to  the  primary  axes  of  the  plot.  This  allows  a  closer  examination  of 
the  higher  order  effects.  The  second-order  contribution  to  SM  comes  from  the  following 
tenn  of  Equation  (1.22): 

l-^(3a0)2r  (2.13) 

Figure  5  shows  the  non-rotated  covariance  ellipsoids  propagated  to  two  orbits 
using  the  same  standard  deviations  given  in  (2.11).  The  plot  on  right  of  Figure  5  is 
zoomed  in  to  show  the  second-order  ellipse  curvature. 
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Non-rotated  Ellipsoid  After  2  Orbits  Non-rotated  Ellipsoid  After  2  Orbits 


Figure  5.  Non-Rotated  Covariance  Ellipsoids  after  Two  Orbits 


The  higher  order  contribution  to  SM  is  small  compared  to  the  total  value  of  8M 
but  it  is  large  compared  to  the  ellipsoid  minor  axis.  A  property  of  linearly  propagated 
covariance  is  that  the  area  of  the  linear  ellipsoid  is  constant  with  time.  The  covariance 
does  not  grow  in  the  5a  direction,  so  as  it  grows  in  the  8M  direction,  the  ellipse  becomes 
thinner  with  time  along  the  minor  axis. 

To  examine  the  higher  order  effects  it  is  necessary  to  determine  a  coordinate 
system  where  the  first-order  covariance  is  diagonal.  This  will  rotate  the  axes  of  the 
linearly  propagated  covariance  to  the  primary  axes  of  the  plot.  To  do  so,  we  use  the 
following: 

Jy  =  RJx  (2.14) 


where  8x  is  the  state  given  Equation  (1.22)  and  5 y  is  the  rotated  state: 


8y  = 


ySu 

_ySa  _ 


(2.15) 


R  is  the  clockwise  rotation  matrix: 
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R  = 


cos  # 
-sin# 


sin# 

cos# 


The  covariance  rotation  is  accomplished  by: 

P  =  RPR 

where  PY  is  the  linearly  propagated  covariance 


and  the  rotated  covariance  P ,  is  given  by: 


Pv  = 


p  p 

rvU  rx 


y  12 


v  12 


1  y  22 


(2.16) 


(2.17) 


(2.18) 


(2.19) 


which  gives: 

Pyu=\(P,u+P,v)  +  \(P,»-P,22y°s(20)  +  P,nM2<>) 

P,n=\(P,u  +P^)-\(P,u -P,n)cos(2g)-PxUsm(2g)  (2.20) 

P,,  =  P,,cos(2<?)-  t(P,  P  .. ) sin i  20) 


To  obtain  the  coordinate  system  where  the  covariance  is  diagonal,  we  set  P  n  =  0 
and  solve  for  the  rotation  angle  # : 


#  =  —  arctan 
2 


2  P 


x\2 


P  -  P 

V  1  xl  1  1  x22  J 


(2.21) 


or 
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(2.22) 


0  =  —  arctan 
2 


2( 

K* 

)“-Km)2 

The  rotation  angle  obtained  from  Equation  (2.22)  is  then  used  with  the  rotation  matrix 
given  by  Equation  (2.16). 

C.  PROBABILITY  DISTRIBUTION 

1.  Overview 


This  section  examines  the  cumulative  probability  distribution  using  the  first  and 
second-order  approximation  as  well  as  the  exact  solution  for  covariance  propagation. 
The  results  of  the  probability  distribution  are  compared  to  the  theoretical  cumulative 
distribution  function  to  assess  the  perfonnance  of  the  approximation.  The  Gaussian 
multivariate  probability  density  function  is: 


/« 


(2.23) 


where  N  is  the  dimension  of  the  state.  The  probability  density  function  as  a  function  of 
number  of  standard  deviations,  k  is: 


m = 


i 

2f,r(U 


(2.24) 


The  Gamma  function  T(x)  is  given  by: 

Y{x)  =  ^ux'e  Udu 

The  probability  distribution  function  as  a  function  of  k  is  then  given  by: 


(2.25) 


21 


F(k)  = 


1 


(2.26) 


— i 

2 2  r 


o 


For  our  purposes  N=2.  The  theoretical  probability  density  function  and  probability 
distribution  function  for  a  two  dimensional  system  are: 


-P 

f{k)  =  he  2 


-P 

F(k)  =  l-e  2 


(2.27) 


2.  Probability  Distribution  Results 


For  the  probability  distribution,  five  evenly  spaced  points  were  investigated  over  a 
range  of  uncertainty  of  the  normalized  semi-major  axis.  The  range  used  for  this  analysis 
is  lxlCT4  to  3xl0'3  corresponding  to  a  range  of  semi-major  axis  uncertainty  of  0.7  km  to 
21  km.  This  is  summarized  in  Table  3.  Initially  the  uncertainty  in  mean  anomaly  is  left 
at  0.01°  through  the  range  of  semi-major  axis  uncertainty. 


Range  of  Uncertainty  in  Normalized  Semi-Major  Axis 


<7. s'. 


1.000x10" 


1.250x10" 


1.550x10" 


2.275x10' 


3.000x10' 


Table  3.  Normalized  Semi-Major  Axis 


Given  Equation  (2.27),  at  k=  3  for  the  theoretical  linear  system,  F(k)  is  0.9889. 
That  is,  a  2D  system  at  3cr  standard  deviation  should  theoretically  capture  98.9%  of  the 
distribution.  In  Table  4,  the  data  presented  shows  the  value  of  k  when  98.9%  of  the 
distribution  is  captured  for  each  of  the  methods  of  propagation.  If  the  approximation  is 
adequately  capturing  the  dynamics,  the  k  @  98.9%  row  should  remain  close  to  3.  If  it 
does  not,  this  indicates  more  standard  deviations  are  required  to  capture  the  system 
dynamics.  The  values  of  F(k)  at  k=  3  for  the  first  and  second-order  approximation  as  well 
as  the  exact  propagation  are  is  also  shown  in  Table  4.  This  allows  us  to  evaluate  how 
much  the  propagation  degrades  with  time. 
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During  this  analysis,  it  was  discovered  that  the  o  .  =1.0x1  CP  /  a5M  =0.0  Tease 
showed  excellent  agreement  in  the  first-order  approximation.  This  is  due  to  the  higher 
order  effects  being  a  function  of  5a  and  the  percentage  change  between  the  first  and 
second-order  approximations  being  smaller  with  larger  aSM.  The  results  presented  in  this 
section  start  at  crf;,  =  8. 25x1 0~4  or  5.775  km.  The  behavior  corresponding  to  the 
<7ga,  =  1 .0x1 0  4  case  will  be  presented  later  along  with  discussion  on  the  relative  behavior 
between  <r .  and  cr,,,. 

Sa  OM 


Probability  Distribution  -  trSa-=  8.25xl0'4 

Orbits 

2 

4 

6 

8 

10 

20 

40 

60 

80 

100 

Linear 

k  @  98.9% 

3.05 

3.21 

3.51 

3.80 

4.19 

6.51 

11.96 

17.56 

23.30 

29.10 

F(k)  @  k= 3 

0.987 

0.983 

0.974 

0.964 

0.951 

0.882 

0.770 

0.684 

0.619 

0.569 

2nd 

Order 

k  @  98.9% 

3.04 

3.16 

3.37 

3.54 

3.73 

4.45 

4.91 

5.00 

5.05 

5.07 

F(k)  @  k= 3 

0.988 

0.985 

0.979 

0.973 

0.967 

0.949 

0.940 

0.937 

0.936 

0.935 

Exact 

k  @  98.9% 

3.06 

3.19 

3.40 

3.59 

3.79 

4.54 

4.95 

5.06 

5.11 

5.13 

F(k)  @  k= 3 

0.987 

0.984 

0.977 

0.972 

0.966 

0.947 

0.938 

0.934 

0.933 

0.933 

Table  4.  Probability  Distribution  -  crSa  =  5.775  km 


One  can  see  in  Table  4  that  the  first-order  approximation  begins  to  significantly 
degrade  by  20  orbits.  This  is  further  illustrated  in  Figure  6.  The  dashed  red  line 
represents  the  theoretical  probability  distribution  from  Equation  (2.27).  The  remaining 
curves  are  the  probability  distribution  plotted  at  20  orbit  increments  out  to  100  orbits. 
The  linear  propagation  curves  flatten  out  with  propagation  time.  This  indicates  a  larger 
standard  deviation  is  needed  to  capture  the  distribution. 
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Figure  6.  Probability  Distribution  Comparison,  100  Orbits  -  oSa  =  5.775  km 


Tables  5  through  7  are  the  summary  of  the  probability  distribution  for  the 
remainder  of  the  range  of  uncertainty  in  semi-major  axis.  The  remaining  probability 
distribution  comparison  plots  are  left  for  the  appendix.  It  can  be  seen  in  Tables  5  through 
7  that  as  we  increase  the  initial  uncertainty  in  semi-major  axis,  leaving  the  initial 
uncertainty  in  mean  anomaly  at  0.01°,  the  accurateness  of  the  first-order  approximation 
decreases. 


Probability  Distribution  -  a6a<  =  1.55xl0'3 

Orbits 

2 

4 

6 

8 

10 

20 

40 

60 

80 

100 

Linear 

k  @  98.9% 

3.68 

5.11 

6.85 

8.78 

10.71 

20.59 

40.94 

61.40 

81.69 

102.06 

F(k)  @  k= 3 

0.969 

0.924 

0.875 

0.830 

0.792 

0.646 

0.497 

0.417 

0.367 

0.331 

2nd 

Order 

k  @  98.9% 

3.44 

4.09 

4.51 

4.73 

4.85 

5.02 

5.09 

5.10 

5.10 

5.10 

F(k)  @  k=  3 

0.976 

0.959 

0.948 

0.943 

0.941 

0.936 

0.935 

0.935 

0.935 

0.935 

Exact 

k  @  98.9% 

3.49 

4.18 

4.59 

4.79 

4.93 

5.08 

5.13 

5.15 

5.13 

5.12 

F(k)  @  k=  3 

0.974 

0.956 

0.946 

0.941 

0.939 

0.933 

0.932 

0.932 

0.933 

0.933 

Table  5.  Probability  Distribution  -  cr Sa  =  10.85  km 

24 


Probability  Distribution  -  oSa.  =  2.275xl0"3 

Orbits 

2 

4 

6 

8 

10 

20 

40 

60 

80 

100 

Linear 

k  @  98.9% 

5.37 

9.38 

13.61 

17.82 

22.15 

44.14 

87.93 

131.9 

175.6 

219.1 

F(k)  @  k=  3 

0.916 

0.818 

0.740 

0.682 

0.630 

0.482 

0.355 

0.292 

0.254 

0.226 

2nd 

Order 

k@  98.9% 

4.16 

4.77 

4.92 

4.99 

5.03 

5.10 

5.10 

5.11 

5.10 

5.09 

F(k)  @  k=  3 

0.957 

0.942 

0.939 

0.937 

0.936 

0.935 

0.935 

0.935 

0.935 

0.935 

Exact 

k  @  98.9% 

4.26 

4.85 

5.00 

5.06 

5.10 

5.14 

5.13 

5.13 

5.14 

5.13 

F(k)  @  k=  3 

0.954 

0.940 

0.936 

0.934 

0.933 

0.932 

0.932 

0.933 

0.933 

0.933 

Table  6.  Probability  Distribution  -crSa  =  15.925  km 


Table  7  shows  that  with  an  initial  uncertainty  in  the  semi-major  axis  equal  to 
3x10'  after  2  orbits  F(k)  at  k= 3  has  already  degraded  to  84%  indicating  that  the  linear 
approximation  is  no  longer  representing  a  significant  portion  of  the  system  dynamics. 
Figure  7  also  illustrates  this.  In  addition  to  the  probability  distribution  plot,  Figure  7 
shows  the  uncertainty  ellipsoid  propagation  for  comparison. 


Probability  Distribution  -  a5a,  =  3xl0'3 

Orbits 

2 

4 

6 

8 

10 

20 

40 

60 

80 

100 

Linear 

k  @  98.9% 

8.30 

15.57 

23.07 

30.61 

38.43 

76.44 

152.6 

228.4 

304.2 

379.9 

F(k)  @  k= 3 

0.842 

0.709 

0.621 

0.556 

0.507 

0.380 

0.271 

0.222 

0.191 

0.173 

2nd 

Order 

k  @  98.9% 

4.70 

4.97 

5.03 

5.07 

5.10 

5.10 

5.09 

5.08 

5.08 

5.08 

F(k)  @  k=  3 

0.944 

0.938 

0.936 

0.935 

0.935 

0.935 

0.935 

0.935 

0.935 

0.935 

Exact 

k  @  98.9% 

4.74 

5.03 

5.10 

5.14 

5.14 

5.14 

5.13 

5.13 

5.13 

5.13 

F(k)  @  k= 3 

0.941 

0.935 

0.933 

0.933 

0.932 

0.933 

0.933 

0.933 

0.933 

0.933 

Table  7.  Probability  Distribution  -  <rSa  =21  km 


It  is  noteworthy  to  see  that  throughout  the  entire  range  of  uncertainty  in  semi¬ 
major  axis,  the  second-order  approximation  is  consistent  out  to  100  orbits  and  closely 
matches  the  exact  probability  distribution.  This  indicates  that  for  the  two  dimensional 
semi-major  axis/mean  anomaly  state,  the  second-order  approximation  is  very  close  to  the 
true  system  dynamics. 
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Figure  7.  Probability  Distribution  Comparison,  10  Orbits  and  Gaussian  Distribution 

Propagated  for  2  Orbits  -  <JSa  =  21  km 


3.  Summary  of  Probability  Distribution 

The  probability  distribution  is  an  excellent  way  to  examine  how  well  a  model  is 
fully  capturing  the  dynamics  of  a  system.  The  results  indicate  that  with  larger  initial 
uncertainty  in  semi-major  axis,  the  first-order  approximation  fails  much  quicker  than  with 
smaller  initial  uncertainty. 

The  results  for  the  o\;.  =1.0x10  4  were  not  presented  in  this  section  and  are 
included  later  in  a  section  that  evaluates  the  relative  behavior  of  <jga ,  and  crm. 
Additionally,  only  variations  in  <jga,  were  presented.  The  results  discussing  the  relative 
behavior  will  show  that  with  a  larger  initial  crSM  relative  to  a  .,  there  is  less  of  a 
percentage  difference  between  the  first-order  and  higher  order  solutions. 
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IV.  PERFORMANCE  METRICS 


A.  INTRODUCTION 

There  has  been  significant  work  done  in  the  area  of  nonlinear  perfonnance 
metrics.  This  chapter  is  an  examination  of  previous  work  in  the  area  of  performance 
metrics  of  prediction  error  and  uncertainty  effects  as  well  as  an  examination  of  the 
relative  behavior  between  initial  uncertainty  in  mean  anomaly  and  semi-major  axis.  The 
ultimate  goal  is  to  detennine  if  a  metric  can  be  developed  to  predict  when  ignored 
nonlinear  effects  on  covariance  propagation  have  rendered  an  approximation  no  longer 
useful. 

B.  PREDICTION  ERROR  COST  FUNCTION 
1.  Overview 

Horwood,  et  at.  [11]  developed  a  metric  for  uncertainty  consistency  called  the 
prediction  error  to  help  identify  errors  in  nonlinear  and  linear  propagation  techniques. 
The  prediction  error  is  given  by: 

PS{t)  =  J" /?j  (*,t)p2  (x,t)r/x  (3.1) 

where  px  and  p2  are  the  probability  density  functions  for  two  satellite  states  evaluated  at  a 
common  time  integrated  over  the  entire  state  space.  The  cost  metric  is  defined  by: 

c(t)  =  -In  V£ (?)  (3.2) 


or 

c(0  =  ^(xi  -x2)r  (pi  P2 ) ~ 1  ( x i  -x2)  +  ^ln[det(2^(P1  +P2))]  (3.3) 

The  first  term  in  Equation  (3.3)  is  one  half  of  the  square  of  the  Mahalanobis  distance  [12] 
and  the  second  tenn  is  the  detenninant  of  the  covariance.  It  will  be  useful  to  analyze 
each  of  these  terms  separately: 


27 


(3.4) 


ci  (0 =^-(xi  -x2)r(pi +p2r  (xi  -x2) 

c2(0=|in[det(2;r(pi+p2)). 


where  the  total  cost,  c(7)  =  cx  (t)  +  c2  (f).  Additionally,  for  purposes  of  this  thesis, 
although  unrealistic,  we  will  be  using  a  “truth”  orbit  for  orbit  2.  Therefore,  the  state 
deviation  x2  and  uncertainty  P2  are  both  zero.  Equation  (3.4)  reduces  to: 


ci(0=|xfpr'xi 

:  (0  =  7T ln  [det  (2^rP[ )] 


(3.5) 


with  the  sum  of  cx  and  c2  being  the  entire  cost  function.  The  first  tenn  of  the  cost  function 
uses  the  Mahalanobis  distance,  so  prior  to  evaluating  the  cost  function,  we  will  examine 
the  MD  with  respect  to  the  state. 

2.  Mahalanobis  Distance 

The  Mahalanobis  distance  is  a  way  to  determine  the  statistically  weighted 
difference  between  two  points.  The  Mahalanobis  distance  is  given  by: 

k2  =  SxTP~lSx, 

where  (3.6) 

8x  =  x(t)-x(t) 


where  the  “bar”  over  the  state  indicates  a  value  on  the  reference  orbit.  Recall  from 
Equation  (1.13)  that  Sa  =  Sa0  and  our  reference  orbit  is  the  truth  orbit.  Therefore: 


a-a 

Sa0 

8M 

(3.7) 
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1st  Order  MD  for  100  Orbits 
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Figure  8.  Mahalanobis  Distance  -  First  and  Second-Order  Approximations 


The  Mahalanobis  distance  for  the  first  and  second-order  propagation  is  plotted  in 
Figure  8.  For  a  linear  system,  the  Mahalanobis  distance  is  constant.  Our  system  is 
nonlinear  but  uses  approximations  to  propagate  the  uncertainty.  As  one  would  expect  for 
the  initial  conditions  SaQ  -  3<r  .  and  SM {]  -  0  the  MD  starts  at  3  and  grows  larger  with 
time.  For  the  linear  propagation,  the  MD  starts  at  3  and  grows  infinitely  large  with  time. 
Interestingly,  for  the  second-order  approximation,  the  MD  approaches  an  asymptote. 
Flow  quickly  the  second-order  MD  approaches  the  asymptotic  value  is  proportional  to 
<7  ..  Smaller  cr  .  results  in  smaller  change  in  MD  suggesting  the  first  and  second-order 
approximations  show  good  agreement  with  the  exact  propagation.  It  will  be  shown  later 
that  this  is  also  a  function  of  the  relative  behavior  between  initial  uncertainty  in  Sa*0  and 
SM0.  Overall,  the  departure  of  the  Mahalanobis  distance  from  the  initial  value  at  initial 
time  is  a  method  of  indicating  how  well  your  approximation  is  capturing  the  full 
nonlinear  dynamics. 
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3. 


Prediction  Error  Cost  Function  Results 


The  cost  function  will  be  used  to  compare  the  results  of  the  first  and  second-order 
approximations.  In  order  to  obtain  the  initial  covariance,  we  must  define  our  initial 
uncertainty.  The  uncertainty  in  mean  anomaly  is  left  at  0.01°.  Equation  (1.26)  shows 
that  it  is  the  uncertainty  in  the  semi-major  axis  or  the  radial  direction  that  causes  the  mean 
anomaly  error  or  in-track  error  to  grow  with  time.  Due  to  this,  we  will  consider  a  range 
of  uncertainties  in  semi-major  axis.  To  get  an  accurate  comparison  with  previous  results, 
the  non-dimensional  semi-major  axis  uncertainties  considered  are  lxl O'4  to  3x10° 
corresponding  to  a  range  of  uncertainty  from  0.7  km  to  21  km.  One  can  see  from 
Equation  (3.3)  that  the  initial  cost  at  time  zero  is  dependent  on  the  selection  of  the  initial 
conditions.  The  initial  condition  used  for  prediction  error  cost  analysis  is  Sa0  =  3o\;, 
and  SM0  =0.  A  different  initial  condition  would  result  in  a  different  value  for  initial  and 
final  cost  and  ultimately  a  different  total  cost.  The  choice  of  initial  condition  is  based 
Sa*0  =  3o\a.  being  the  “worst  case”  from  a  range  of  uncertainty  bounded  by  the  3cr 
ellipse.  The  total  cost  for  10  orbits  is  plotted  in  Figure  9. 


Total  1st  OrSer  Cost  Function  for  10  Orbrts 


Figure  9.  Total  Cost  Metric  for  10  Orbits  -  First  and  Second-Order  Approximation 
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The  results  of  the  total  cost  as  well  as  the  Mahalanobis  distance  for  the  first-order 
approximations  are  summarized  in  Table  8.  Propagation  times  are  set  at  10  as  well  as 
100  orbits  for  comparison.  Interestingly,  for  a  small  initial  uncertainty,  or  a ga.  =  1x10  4, 
there  is  very  little  change  in  total  cost  or  MD  over  both  the  10  and  100  orbit  propagation 
times.  This  initially  would  indicate  that  the  uncertainty  consistency  is  well  represented 
for  smaller  Like  the  Mahalanobis  distance,  the  performance  of  the  <Jga.  =  1x10  4  case 
is  in  reality  maximizing  perfonnance  of  the  relative  behavior  between  <JSM  and  a  .. 

The  c2  term  of  the  first-order  approximation  requires  some  discussion.  The  second 
tenn  of  the  cost  is  a  function  of  the  determinant  of  the  covariance.  Alfriend  [1]  proved 
that  for  Hamiltonian  systems,  the  volume  of  the  uncertainty  ellipsoid  is  constant  with 
time.  In  other  words,  the  determinant  of  the  covariance  is  constant.  This  was  also 
demonstrated  during  the  covariance  ellipsoid  propagation  with  the  “thinning”  of  the 
ellipsoid  as  it  grew  in  the  5M  direction.  Due  to  the  constant  area  of  the  covariance 
ellipsoid  propagated  under  the  first-order  approximation,  the  difference  in  the  c2  term  or 
the  total  cost  in  c\  over  the  propagation  time  is  zero. 


1st  Order  Approximation 

10  Orbits 

Total  Cost 

Term  1  Cost 

Term  2  Cost 

MD 

1.000x1  O'4 

1.8440x1  O’3 

1.8440xl0'3 

0 

3.000615 

8.250x1 0"4 

8.4991 

8.4991 

0 

5.098846 

1.550xl0‘3 

105.3632 

105.3632 

0 

14.823173 

2.275  xlO’3 

486.5164 

486.5164 

0 

31.337402 

3. 000x1  O’3 

1463.7658 

1463.7658 

0 

54.189773 

100  Orbits 

Total  Cost 

Term  1  Cost 

Term  2  Cost 

MD 

1.000x1  O'4 

0.1844 

0.1844 

0 

3.060849 

8.250xl0"4 

849.9114 

849.9114 

0 

41.337911 

1.550xl0"3 

1.0536xl04 

1.0536xl04 

0 

145.195203 

2.275xl0‘3 

4.8652xl04 

4.8652xl04 

0 

311.949164 

3. 000x1  O’3 

1.4638x10s 

1.4638xl05 

0 

541.075005 

Table  8.  First-Order  Approximation  Cost  and  Mahalanobis  Distance 
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The  results  of  the  total  cost  for  the  second-order  approximation  are  summarized  in 
Table  9.  As  one  would  expect,  overall  the  total  cost  for  the  first-order  approximation  is 
greater  than  the  second-order  approximation  indicating  greater  uncertainty  consistency 
with  a  second-order  approximation.  However,  this  does  not  hold  with  the  <Jga.  =  lxl  0  4 
case.  The  total  cost  for  this  case  is  very  close  for  the  first  and  second-order 
approximation.  The  case  with  initial  uncertainties  of  a  Sa,  =  1x1 0  4  and  aSM  =  0.01°  show 
good  agreement  between  the  first  and  second-order  approximations.  This  agrees  with  the 
results  obtained  from  the  probability  distribution  and  is  also  due  to  the  relative  behavior 
in  initial  uncertainties. 

As  stated  previously,  a  property  of  a  linearly  propagated  covariance  is  that  the 
volume  is  constant  with  time.  The  tenn  c2  is  simply  the  volume  of  the  covariance  scaled 
by  a  constant.  The  condition  of  the  covariance  ellipsoid  volume  being  constant  with  time 
does  not  hold  for  the  second-order  approximation.  One  can  see  there  is  a  contribution  for 
the  Term  2  cost  at  a  .  =  1.0x10  .  This  actually  results  in  a  slightly  higher  cost  for  the 
second-order  approximation  for  this  case. 


2nd  Order  Approximation 

10  Orbits 

C»a 

Total  Cost 

Term  1  Cost 

Term  2  Cost 

MD 

l.OOOxlO"4 

1.8894xl0'3 

1.8438xl0'3 

4.5560xl0'5 

3.000615 

8.250xl0'4 

6.1524 

5.9763 

0.1761 

4.577401 

1.550xl0'3 

17.7490 

16.8319 

0.9171 

6.531751 

2.275  xlO'3 

20.7644 

19.1468 

1.6176 

6.877041 

3.000xl0'3 

21.7236 

19.5661 

2.1575 

6.937740 

100  Orbits 

C»a 

Total  Cost 

Term  1  Cost 

Term  2  Cost 

MD 

l.OOOxlO’4 

0.1873 

0.1827 

0.0045 

3.060305 

8.250xl0’4 

21.5507 

19.6676 

1.8831 

6.952358 

1.550xl0'3 

23.1276 

19.9940 

3.1336 

6.999144 

2.275xl0'3 

23.8234 

19.9231 

3.9003 

6.989003 

3.000xl0'3 

24.2819 

19.8285 

4.4534 

6.975459 

Table  9.  Second-Order  Approximation  Cost  and  Mahalanobis  Distance 


32 


4.  Summary  of  Cost  Function  and  Mahalanobis  Distance 


Both  the  prediction  error  and  the  Mahalanobis  distance  give  comparison  between 
an  initial  and  propagated  state.  The  prediction  error  takes  into  account  higher  order 
statistics  with  the  c2  term.  Neither  metric  provides  definitive  information  or  give  a 
consistent  means  of  determining  when  the  linear  approximation  is  no  longer  valid.  The 
prediction  error  cost  function  is  very  useful  for  determining  perfonnance  gains  when 
using  higher  order  filters,  but  overall  doesn’t  answer  the  problem  posed  in  this  thesis. 

The  Mahalanobis  distance  gives  a  statistical  difference  between  an  initial  point 
and  the  same  point  after  propagation.  Alfriend  [1]  has  shown  that  the  Mahalanobis 
distance  can  be  used  for  uncorrelated  track  association  in  a  more  statistically  valid 
manner.  With  these  methods  however,  the  selection  of  a  specific  MD  as  a  metric  is  a 
balance  between  accurate  correlations  and  false  correlations.  Previous  work  by  Hill  et  al. 
[3]  suggest  a  MD  of  less  than  10  can  be  used  for  UCT  association  with  a  reasonable 
amount  of  confidence.  There  is  no  definitive  total  change  in  MD  that  identifies  when  the 
first-order  approximation  is  no  longer  valid. 

C.  UNCERTAINTY  METRICS 

In  addition  to  the  work  with  prediction  error  cost  functions,  there  have  been  other 
studies  into  determining  metrics  to  assess  and  compare  performance  of  higher  order 
propagation  techniques.  In  an  attempt  to  detennine  if  these  metrics  will  reveal  any 
information  into  the  validity  of  a  first-order  approximation,  we  will  examine  these 
metrics  further. 

Park  and  Scheeres  developed  a  metric  to  assess  the  perfonnance  of  the  State 
Transition  Tensor  method  of  mapping  uncertainties.  It  was  found  that  when  applied  to 
the  state  used  in  this  thesis,  this  metric  displayed  some  singularities  close  to  the  initial 
time,  r  =  0.  Due  to  this,  it  was  required  to  slightly  modify  the  Park  and  Scheeres  metric. 
The  results  of  the  modified  metric  are  presented  in  this  section.  To  begin,  a  review  of  the 
Park  and  Scheeres  metric  is  conducted. 
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1. 


Park  and  Scheeres  Nonlinearity  Metric 


The  Park  and  Scheeres  [7]  local  nonlinearity  index  was  proposed  to  assist  with 
determining  the  sufficient  order  of  higher-order  solutions.  This  local  nonlinearity  index 
is  defined  as: 
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)-Sx* 
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Sx*( 

t;Sx°k,t( 
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(3.8) 


where  Sx is  the  mth  order  approximation  of  Equation  (1.8),  Sx°k  is  the  Mi 
initial  condition  for  Sx  chosen  from  the  boundary  of  the  initial  confidence  region,  i.e.  on 
the  surface  of  the  ka  ellipse,  and  Sx*  is  the  exact  solution  of  Equation  (1.8).  Einstein’s 
summation  notation  is  used  and  n  is  the  dimension  of  the  state  while  N  is  the  number  of 
initial  conditions. 
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Figure  10.  Park  and  Scheeres  Nonlinearity  Index 


Upon  examination  of  this  metric,  it  was  discovered  that  for  the  state  used  in  this 
thesis,  there  were  singularities  due  to  the  denominator  of  Equation  (3.8)  occasionally 
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being  very  close  to  zero,  particularly  near  time  r  =  0 .  The  singularities  are  illustrated  in 
Figure  10.  The  results  in  Figure  10  were  obtained  using  Equation  (3.8)  with  1000  initial 
conditions  sampled  on  the  3cr  ellipsoid  with  ( r  .  =  0.003  and  <JSM  =  0.01°. 


The  steady  state  value  of  rj  for  this  particular  set  of  initial  conditions  and  initial 
uncertainty  is  actually  1.1255x10'  for  the  first-order  approximation  and  1.1823x10'  for 
the  second-order  approximation.  One  can  see  the  singularities  at  the  beginning  of  the 
propagation  are  distorting  the  curve.  The  following  will  derive  the  limit  of  the  index  in 
order  to  determine  what  the  steady  state  behavior  is. 

Recall  that  the  expression  for  the  exact  solution  of  8M  is: 


*  Vi  *\  —3/2 

8M  =2tt  (a  )  -1  r  +  8Mt 


/  *  \  —  3/2 

The  Taylor  series  expansion  of  (a  j  is  given  by: 


j= i 

1  d(ajm 

Cj  j\  da 


(3.10) 


therefore  the  exact  solution  for  8M  becomes: 


oo 

SM*  =  IttJc,  (Sa*  )J t  +  SM0 


(3.11) 


and  the  mth  order  approximation  for  the  change  in  mean  anomaly  is: 


m 

8Mm  =  2nY4cj  (Sa)JT  +  8M0 


(3.12) 


The  c to  m=3  are: 


(3.13) 


Although  the  nonlinearity  index  is  a  function  of  the  entire  state,  5a  is  constant. 
Therefore,  only  the  expression  for  5M  need  be  considered  for  the  index.  Substituting 
Equation  (3.11)  and  (3.12)  into  Equation  (3.8),  the  local  nonlinearity  index  becomes: 
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Consolidating  the  sigma  notation  in  the  numerator  gives: 
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With  the  results  in  Equation  (3.15),  there  are  two  observations  about  the 
performance  of  this  metric  applied  to  our  state.  One  is  that  when  5M°k  and 

00 

2/T^cv  [Sak^  r  are  equal  and  of  opposite  sign,  77'"  will  be  infinite.  These  singularities 

7= 1 

often  appear  near  r  =  0.  Due  to  this  it  is  necessary  to  slightly  modify  the  local 

00  j 

nonlinearity  index.  The  second  observation  is  that  eventually  the  X0/  ( 5ak  j  r  term  will 

0  7=1 

dominate  5Mk.  Once  this  occurs,  we  can  express  the  steady  state  value  reached  by 

Equation  (3.15)  as: 
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2.  Proposed  Performance  Index  #1 


In  order  to  mitigate  the  singularity  behavior  with  the  Park  and  Scheeres  local 
nonlinearity  index,  the  problem  with  the  tenn  in  the  denominator  being  zero  needs  to  be 
addressed.  The  first  of  two  proposed  perfonnance  metrics  is  given  by: 
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It  is  important  to  note  that  <7 x  must  be  dimensionless  or  each  term  in  8x  must  have  the 
same  dimension.  Again,  due  to  the  fact  that  8  a  is  constant,  we  can  simplify  the 
expression  for  77  and  express  it  in  terms  of  8M.  Equation  (3.17)  becomes: 
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Equation  (3.18)  will  start  at  zero  and  asymptotically  approach  the  following 

value: 


77 


m 


sup 

k=l,...,N 


(3.19) 


Consider  the  ktj  ,  to  be  the  maximum  value  of  8a.  This  occurs  when  8Mn  =  0.  This 

da 

leaves: 
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where  the  Cj  come  from  the  Taylor  series  expansion  of  I  a  )  about  the  reference  value 
a*  =  1  given  in  Equation  (3.10). 
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1st  Order  Case 
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Figure  11.  Proposed  Performance  Metric  #1-10  Orbits 


Figure  10  shows  Performance  Index  #1  for  the  range  of  <r  ..  The  uncertainty  in 
mean  anomaly  is  left  at  0.01°.  One  can  see  that  for  the  a  .  =  1 .0x1 0  4  case,  the  first-order 
is  still  showing  some  discontinuities  during  initial  propagation  time  and  settling  on  the 
theoretical  value.  It  is  unclear  what  is  causing  these  discontinuities.  Based  on  the  results 
of  this  metric,  it  provides  a  measure  of  increased  performance  with  a  higher  order 
solution,  but  does  not  provide  information  on  when  the  first-order  approximation  is  no 
longer  valid.  Due  to  this,  further  investigation  into  the  discontinuities  is  not  warranted. 
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Table  10.  Proposed  Performance  Metric  #1-10  Orbits 
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3.  Proposed  Performance  Index  #2 


The  second  of  the  proposed  perfonnance  metrics  relies  on  the  initial  uncertainty 
in  mean  anomaly.  It  is  given  by: 
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Using  Equation  (3.11)  and  (3.12),  the  theoretical  behavior  for  Metric  #2  can  be 
expressed  analytically  in  the  following  manner: 
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Given  Equation  (3.24)  and  points  selected  from  the  k<j  ellipsoid,  the  maximum 
if1  occurs  when  8M0-  0  and  |<7a  |  is  a  maximum.  From  Equation  (3.24),  the 
theoretical  77'”  can  be  expressed  analytically  as: 
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The  numerical  results  from  Equation  (3.21)  are  shown  in  Figure  12  for  10  orbits. 
As  with  previous  examples,  the  range  of  ct,.;  is  lxlO'4  to  3xl0'3  and  <rSM  is  0.01°. 
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1st  Order  Case 


2nd  Order  Case 


Figure  12.  Proposed  Performance  Metric  #2-10  Orbits 


Table  11  summarizes  the  results  of  the  numerical  versus  theoretical  comparison 
for  Metric  #2.  As  expected,  the  total  change  in  the  metric  is  larger  for  the  first-order 
approximation.  One  can  see  there  is  good  agreement  between  the  numerical  results 
obtained  through  Equation  (3.21)  and  the  predicted  results  from  Equation  (3.25). 
Interestingly,  the  behavior  of  Metric  2  is  consistent  throughout  the  range  of  uncertainty  in 
5a  and  does  not  appear  to  be  effected  by  the  relative  behavior  of  initial  uncertainties  in 
mean  anomaly  and  semi-major  axis. 
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4.  Summary  of  Metric  Performance 


The  local  nonlinearity  index  developed  by  Park  and  Scheeres  was  developed  to 
decide  the  sufficient  order  of  a  higher  order  solution.  The  original  metric  as  well  as 
Index  #1  represent  the  percent  difference  between  the  approximation  and  the  true 
solution.  This  is  very  useful  if  one  defines  a  threshold  of  accurateness  for  an 
approximation.  Ultimately,  the  metrics  do  not  give  definitive  infonnation  as  to  when  the 
first-order  approximation  is  no  longer  representing  the  system  dynamics. 

D.  RELATIVE  BEHAVIOR  OF  INITIAL  UNCERTAINTIES 

1.  Overview 

Thus  far  it  has  been  suggested  that  for  the  <J Sa.  =0.0001  and  a5M  =0.01°  case, 
there  is  very  good  agreement  over  extended  periods  of  propagation.  This  prompted 
further  investigation  into  the  relationship  between  the  initial  semi-major  axis  uncertainty 
and  the  mean  anomaly. 

This  section  will  examine  the  relative  behavior  between  the  initial  uncertainties 
and  the  higher  order  effects  on  8M .  Additionally,  a  performance  metric  is  presented  that 
predicts  how  long  the  linear  approximation  is  valid  given  a  target  distribution  and  initial 
uncertainties. 

2.  First-order  Approximation  at  Small  Initial  Uncertainty 

Figure  13  shows  the  initial  distribution  of  106  sample  points  propagated  using  the 
complete  system  dynamics,  and  then  the  3cr  uncertainty  ellipse  propagated  using  the  first 
and  second-order  approximations  (green  ellipse  and  blue  ellipse  respectively).  At  100 
orbits,  there  is  very  little  difference  between  the  first  and  second-order  approximations 
and  the  first-order  approximation  still  encapsulates  nearly  99%  of  the  initial  distribution. 
The  initial  3cr  ellipse  for  the  simulation  captured  99.92%  of  the  106  points.  MATLAB’s 
pseudorandom  number  function  randn()  was  used  to  generate  the  initial  distribution.  One 
can  see  that  even  after  100  orbits,  the  linear  approximation  identified  by  the  green  3cr 
ellipse  is  representing  the  true  system  dynamics  fairly  well.  The  blue  ellipse,  which  is 


41 


slightly  distorted,  is  the  second-order  approximation,  but  overall  there  is  very  little 
difference. 


3ct  Covariance  After  100  Orbits 
Distribution  Captured  by  3g  Ellipse  =  0  9892 


y5M 


xIO 


-7 


*  Monte  Carlo  - 1st  Order - 2nd  Order 


Figure  13.  3cr  Uncertainty  Ellipse  and  Gaussian  Distribution  Propagated  for  100  Orbits  - 

=  0.7  km 


Likewise,  Figure  14  shows  the  probability  distribution  curve  using  the  same  initial 
uncertainty  in  semi-major  axis.  At  100  orbits,  there  is  near  perfect  agreement  in  first- 
order,  second-order  and  exact  propagations.  Table  12  summarizes  the  results  for  the 
probability  distribution  out  to  100  orbits. 
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Figure  14.  Probability  Distribution  Comparison,  100  Orbits  -  cts  =  0.7  km 
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Table  12.  Probability  Distribution  -  <JSa  =  0.7  km 


Finally,  looking  at  the  results  of  the  Horwood  Prediction  Error  Cost  Function,  the 
Mahalanobis  distance,  and  the  performance  indices  #1  and  #2,  all  of  these  metrics 
suggest  that  the  &Sa.  =0.0001  case,  there  is  excellent  agreement  between  first-order 
approximations. 

Initially,  given  that  the  equations  of  motion  are  a  function  of  semi-major  axis,  it 


seems  intuitive  that  with  smaller  initial  uncertainty  in  semi -major  axis,  the  first-order 
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approximation  would  remain  valid  for  longer  propagation  time.  However,  the  results  at 
the  <j  .  =0.0001  /  <jsm  =0.001°  case  show  the  first-order  approximation  captures  the 
system  dynamics  out  to  extremely  extended  propagation  times,  past  even  1000  orbits. 
Interestingly  it  was  found  that  if  the  initial  mean  anomaly  uncertainty  was  varied  so  that 
it  was  much  smaller  than  the  initial  uncertainty  in  semi -major  axis,  the  validity  of  the 
approximation  was  not  as  good  at  extended  orbital  periods. 

Figure  15  shows  the  probability  distribution  as  well  as  the  covariance  propagation 
out  to  100  orbits  using  <jga.  =0.0001.  The  uncertainty  in  mean  anomaly  however  is 
changed  from  the  original  0.01°  to  crSM  =  0.001°.  There  is  still  fairly  good  agreement  at 
100  orbits,  but  not  as  good  as  the  asu  =0.01°  case.  The  right  of  Figure  15  shows 
approximately  92%  of  the  distribution  captured  by  the  first-order  approximation  at  100 
orbits  versus  98.9%  at  100  orbits  for  crSM  =0.01°.  This  suggests  that  the  validity  of  the 
first-order  approximation  is  a  function  of  the  initial  mean  anomaly  uncertainty,  which 
was  unexpected. 
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Figure  15.  3cr  Uncertainty  Ellipse  and  Gaussian  Distribution  Propagated  for  100  Orbits  - 

(7Sa=  0.7  km  and  crSM  =  0.001° 
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3. 


First-order  Approximation  Error 


In  this  section,  a  comparison  of  the  first-order  approximation  and  the  exact 
solution  is  made  in  order  to  obtain  an  expression  for  the  validity  of  the  linear 
approximation  as  a  function  of  time  and  initial  uncertainties.  In  doing  so,  this  will  also 
investigate  the  relative  behavior  between  initial  uncertainty  in  mean  anomaly  and  semi¬ 
major  axis  and  the  influence  on  validity  of  propagation  time. 

The  exact  solution  for  difference  in  mean  anomaly  expressed  in  non-dimensional 
tenns  is: 
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The  first-order  approximation  is: 
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Define  y  as  the  first-order  approximation  error: 
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The  difference,  y  is  given  by  Equation  (3.29): 
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Simplified,  Equation  (3.29)  becomes: 
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It  is  notable  that  Equation  (3.30)  is  not  a  function  of  SM0.  This  is  inconsistent 
with  previous  results.  In  particular,  Figures  13  through  15  suggest  the  initial  uncertainty 
in  SM0  can  have  a  significant  influence  on  the  first-order  to  exact  solution  error. 

Figure  16  shows  y  at  1000  orbits  for  a  Sa  range  of  ±60  km  (60  km  representing 
the  approximate  3cr  upper  bound  for  Sa  presented  thus  far  in  the  thesis).  A  significant 
observation  from  Equation  (3.30)  and  Figure  16  is  that  when  Sa  =  Otheny  =  0.  The  fact 
that  the  first-order  approximation  error  is  proportional  to  Sa  is  a  contributing  factor  to 
why  the  results  thus  far  show  good  agreement  for  the  Sa *  =  0.0001.  This  is  only  part  of 
the  solution,  as  Equation  (3.30)  and  Figure  16  fail  to  show  any  correlation  between  initial 
uncertainty  in  mean  anomaly  and  the  first-order  approximation  validity. 
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Figure  16.  First-Order  Approximation  Error 

4.  Linear  Approximation  Performance  Metric 

As  stated  previously,  the  ultimate  goal  is  to  determine  if  there  is  a  performance 
metric  that  can  be  used  to  identify  when  the  nonlinearities  in  the  equations  of  motion 


First-Order  Approximation  Error 
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render  the  first-order  approximation  invalid.  The  following  is  an  attempt  to  analytically 
relate  the  validity  of  the  first-order  approximation  with  initial  uncertainty  in  SM0. 

Given  that  previous  results  show  the  validity  of  the  linear  approximation  is  a 
function  of  aSM ,  the  ka  uncertainties  were  substituted  into  Equation  (3.30)  and  y  was  set 
equal  to  kaSMo  in  order  to  determine  if  this  was  the  time  when  the  linear  approximation 
was  no  longer  capturing  the  system  dynamics. 


\(1+(*<v))’ 


+  —  ka  ,  -1 

2  Sa 


2m  =  ( k(JdMo ) 


(3.31) 


Equation  (3.31)  allows  us  to  obtain  a  time  in  terms  of  number  of  orbits  when  the 
difference  between  the  first-order  approximation  and  exact  solution  is  equal  to  kaSM  . 
This  is  given  by  Equation  (3.32) 


T 


^C(TSMq 


1  1 

(1+(K.)) 

2  n 


(3.32) 


Table  13  summarizes  the  results  for  k=2  using  <rSM  =  0.0  T.  The  same  range  of 
initial  uncertainties  in  semi-major  axis  used  throughout  this  thesis  is  presented  here. 
Additionally  for  reference,  the  probability  distribution  F(k)  at  r  =  r3(J  is  given  as  well  as 
the  value  of  y  at  1000  orbits. 


1.000x1  O'4 

8.250X10-4 

1.550xl0'3 

2.275xl0'3 

3.000xl0'3 

r3-«0(OrbitS) 

494.0 

7.276 

2.067 

0.9617 

0.5545 

F(k)  at  r  =  r3 

SM0 

0.9681 

0.9681 

0.9681 

0.9681 

0.9681 

y  at  r=1000  orbits 

1.0599xl0'3 

7. 195  8x1  O'2 

0.25336 

0.54443 

0.94434 

Table  13.  Summary  of  Performance  of  r3tT 
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As  shown  in  Table  13,  The  first-order  approximation  F(k)  is  still  at  approximately 
96.8%  of  the  distribution  when  3 crSM  is  equal  to  ;/  ,  which  is  still  significant  agreement. 
It  follows  that  setting  kaSM ■  equal  to  y  does  not  give  the  correct  metric  for  assessing  the 
validity  of  the  linear  approximation.  However,  previous  results  show  the  validity  of  the 
linear  approximation  is  a  function  of  kagMa,  so  there  still  is  value  in  examining  Equation 
(3.32). 

To  further  illustrate  the  behavior  of  r3tJ m  ,  Figure  17  shows  Equation  (3.32) 
plotted  over  a  range  of  ±20  km  in  <JSa.  The  behavior  of  r3(T  _  at  crSa  -  0  is  infinite; 
however  the  y-axis  is  scaled  to  1500  orbits  in  order  to  show  the  relevant  behavior.  It’s 
clear  from  Equation  (3.32)  that  with  a  larger  numerator,  the  behavior  of  the  curves 
approaching  infinity  will  extend  outward,  essentially  giving  a  greater  number  of  orbits 
with  larger  values  of  aSa. 
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Figure  17.  r3(^  with  aSUo  =  0.0  F 

With  respect  to  the  results  for  Equation  (3.32)  listed  in  Table  13,  there  is 
interesting  behavior  with  F(k).  One  can  see  that  the  number  of  orbits  obtained  change 


x  versus  q 
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with  different  <JSa,  but  the  probability  distribution,  F(k)  at  r  =  r3cr  ^  does  not  change  with 
different  <7  ..  Using  this  property  will  allow  us  to  detennine  an  expression  for  the 
number  of  orbits  as  a  function  of  a  .  and  crSM  that  target  a  specific  percentage  of  the 
distribution  we  want  the  first-order  approximation  to  capture. 


TF(k)  ~ 


\I(1+(K.-)1' 


+  - 


(KvH 


2  n 


(3.33) 


Given  values  for  a ga.  and  <jSMo  ,  the  number  of  orbits  for  a  target  distribution  can 
be  determined  experimentally  using  either  the  error  ellipse  propagation  or  the  probability 
distribution  calculation.  Equation  (3.33)  can  then  be  solved  for  the  coefficient  x.  We  can 
then  use  Equation  (3.33)  to  determine  the  approximate  number  of  orbits  the  first-order 
approximation  is  good  for  with  a  given  a  target  distribution,  as  a  function  of  a  „  and 
<JSM  .  A  summary  of  the  coefficient,  x  for  various  target  distributions  is  listed  in  Table 
14. 


Target  Distribution 

0.95 

0.9 

0.85 

0.8 

0.75 

0.7 

0.65 

Coefficient 

1.35020 

2.35425 

3.36235 

4.46964 

5.75911 

7.29757 

9.15385 

Table  14.  Summary  of  Coefficients  used  for  Equation  (3.33) 


Ultimately  the  relative  performance  of  the  first-order  approximation  with  respect 
to  oSM  is  due  to  the  percentage  difference  from  the  higher  order  effects.  The  mean 
anomaly  error,  as  well  as  any  higher  order  contributions  to  mean  anomaly  error,  are 
solely  a  function  of  5a.  The  first-order  approximation  for  the  uncertainty  propagation 
has  excellent  agreement  with  the  second-order  and  exact  solution  at  crSa.  =0.0001  and 
osu  =  0.01°  because  the  higher  order  contribution  from  5a  is  very  small  relative  to  <JSM. 


49 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


50 


V.  CONCLUSIONS 


A.  KEY  POINTS  AND  RECOMMENDATIONS 

In  general,  it  was  proven  difficult  to  determine  a  definitive  measure  of  when  the 
covariance  propagated  using  a  linear  approximation  was  no  longer  Gaussian.  The 
cumulative  distribution  comparisons  allow  one  to  assess  the  performance  of  the 
approximation  given  initial  uncertainty.  The  propagation  of  the  covariance  ellipsoid 
along  with  the  Monte  Carlo  simulations  provides  a  good  method  of  visualizing  the 
performance.  Using  these  methods  however,  much  like  using  the  Mahalanobis  distance, 
relies  on  the  user  to  define  a  threshold  that  is  no  longer  acceptable. 

The  metrics  such  as  the  prediction  error  cost  function  and  the  local  nonlinearity 
index  provide  a  way  of  assessing  the  significant  order  of  a  higher  order  solution’s  ability 
to  fully  capture  a  system’s  dynamics.  These  metrics  however  did  not  provide  any 
information  on  which  initial  uncertainties  and  their  size  cause  propagated  covariance  to 
become  non-Gaussian. 

The  performance  of  the  linear  approximation  with  small  initial  uncertainty  in  a  . 
prompted  investigation  into  the  relative  behavior  between  the  initial  uncertainties.  This 
revealed  that  although  the  higher  order  effects  are  a  function  of  8a  ,  they  only  have  a 
significant  impact  when  they  are  large  in  comparison  to  the  uncertainty  ellipse  minor 
axis,  which  is  a  function  of  the  initial  uncertainty  in  mean  anomaly.  During  the 
investigation  into  the  relative  behaviors,  a  linear  approximation  performance  metric  was 
developed.  Although  it  is  only  applicable  to  the  state  as  defined  in  this  thesis,  a  similar 
approach  could  be  used  with  a  more  complex  model. 

B.  AREAS  OF  FURTHER  RESEARCH 

The  system  used  in  this  thesis  is  very  basic,  consisting  of  only  two-body  motion 
in  the  semi-major  axis/mean  anomaly  state  space.  Previous  work  referenced  throughout 
this  thesis  has  shown  that  different  coordinate  systems  can  produce  vastly  different 
results  with  respect  to  nonlinear  behavior  in  first-order  approximations.  Further  research 
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using  different  coordinate  systems  as  well  as  using  more  complex  models  is  warranted. 
Additionally,  the  linear  approximation  perfonnance  metric  developed  in  this  thesis  could 
be  modified  to  be  made  applicable  to  a  real-world  model  and  be  used  to  predict  more 
definitively  how  long  a  linear  approximation  is  valid. 
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APPENDIX  A.  COVARIANCE  EXPANSION  TO  SECOND-ORDER 
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APPENDIX  B.  ADDITIONAL  PLOTS 


Additional  Cost  Function  Plots 


Figure  18.  Cost  Function  Tenn  1  and  2  for  10  Orbits 
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Figure  19.  Cost  Function  Term  1  and  2  for  100  Orbits 
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Figure  20.  Total  Cost  Metric  for  100  Orbits,  First  and  Second-Order  Approximation 


2.  Additional  Probability  Distribution  Plots 
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Figure  21.  Probability  Distribution  Comparison,  10  Orbits  -  crSa  =  0.7  km 
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Figure  22.  Probability  Distribution  Comparison,  10  Orbits  -  <JSa  =  5.775  km 
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Figure  23.  Probability  Distribution  Comparison,  10  Orbits  -  <JSa  =  10.85  km 
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Figure  24.  Probability  Distribution  Comparison,  100  Orbits  -  crSa  =  10.85  km 
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Figure  25.  Probability  Distribution  Comparison,  10  Orbits  -  uSa  =  15.925  km 
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Figure  26.  Probability  Distribution  Comparison,  100  Orbits  -  <JSa  =  15.925  km 


2 

b 


2D 

03 

2D 

O 


Linear  Propagation 

CTga*  =  0003 


k 


2nd  Order  Propagation 


k 


Q 

>, 


XJ 

CD 

n 

P 


0.5 


Exact  Propagation 
cts_*  =  0.003 

Oa 


Xf 


f 


V 


2  4 

k 


—  0  Orbits 

—  2  Orbits 
- 4  Orbits 

6  Orbits 

—  8  Orbits 
10  Orbits 

- Theoretical 


Figure  27.  Probability  Distribution  Comparison,  10  Orbits  -  crSa  =  21  km 
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Figure  28.  Probability  Distribution  Comparison,  100  Orbits  -  crSa  =21  km 


3.  Mean  Anomaly  Error  and  In-Track  Error  Plots 
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- Ojj.  =  0  0001  - Oj,.  *  0  000825  - c^.  =  0  00155  - Ou.  =  0  002275  - o^.  =  0  003 

Figure  29.  Mean  Anomaly  Error,  10  Orbits  -  Sa0  =  3cr . .,  SM0  =  0 
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1st  Order  Mean  Anomaly  Error  for  100  Orbits 


1st  Order  In-Track  Error  for  10  Orbits 


2nd  Order  In-Track  Error  for  10  Orbits 


—  a -  0  0001  - Oj,.  =  0  000825  - =  0  00155  a^.  =  0  002275  - a^.  =  0  003 

Figure  3 1 .  In-Track  Error,  10  Orbits  -  Scin  =  3cr . ,,  SM0  -  0 
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1st  Order  In-Track  Error  for  100  Orbits 


- Oj,.  =  0  0001  - Oj,.  =  0  000825  - c^,  =  0  00155  - au.  =  0  002275  - o^.  =  0  003 

Figure  32.  In-Track  Error,  100  Orbits  -  Sa0  =  3cr ;  SM0  =  0 
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